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ABSTRACT 

Planar flames are intrinsically unstable in open domains due to the thermal 
expansion across the burning front — the Landau-Darrieus instability. This in- 
stability leads to wrinkling and growth of the flame surface, and corresponding 
acceleration of the flame, until it is stabilized by cusp formation. We look at the 
Landau-Darrieus instability for C/0 thermonuclear flames at conditions relevant 
to the late stages of a Type la supernova explosion. Two-dimensional direct 
numerical simulations of both single-mode and multi-mode perturbations using 
a low Mach number hydrodynamics code are presented. We show the effect of 
the instability on the flame speed as a function of both the density and domain 
size, demonstrate the existence of the small scale cutoff to the growth of the in- 
stability, and look for the proposed breakdown of the non-linear stabilization at 
low densities. The effects of curvature on the flame as quantified through mea- 
surements of the growth rate and computation of the corresponding Markstein 
number. While accelerations of a few percent are observed, they are too small to 
have any direct outcome on the supernova explosion. 

Subject headings: supernovae: general — white dwarfs — hydrodynamics - 
nuclear reactions, nucleosynthesis, abundances — conduction — methods: nu- 
merical 



1. INTRODUCTION 

A carbon/oxygen flame propagating outward from the center of a white dwarf is sub- 
jected to a number of instabilities that wrinkle and accelerate the flame. If the flame can 
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accelerate to a significant fraction of the speed of sound, a deflagration alone can account for 
the explosion (see Hillebrandt & Niemeyer 2000 and references therein; Reinecke et al. 2002; 
Gamezo et al. 2003). One such instability, the Landau-Darrieus (LD) instability (Darrieus 
1938; Landau 1944) affects a planar flame front, even in the absence of gravity, driven by the 
thermal expansion across the flame. The LD unstable flame will wrinkle, eventually forming 
cusps in the nonlinear regime that stabilize the flame. The growth rate for a LD unstable 
flame, including the effects of the finite thickness of the flame, //, is 
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(Zeldovich et al. 1985), where a = Pf uc i/p a sh is the density ratio across the flame, k is the 
wavenumber, Ui is the laminar flame speed, and Ma is the Markstein number. The Markstein 
number is a measure of the response of the flame speed, U, to curvature, 

U = + Ma l f ^j , (2) 

where yf is the position of the flame interface and x is the transverse coordinate. We note 
that our sign convention for Ma is opposite that used in Zeldovich et al. (1985) (see Dursi 
et al. 2003 for a discussion of the Markstein number for astrophysical flames). This growth 
rate has been confirmed experimentally (Clanet & Searby 1998) for chemical flames. In this 
paper, we present two-dimensional direct numerical simulations (DNS) of the LD instability 
for low-density C/O flames in white dwarf interiors. We start by briefly reviewing previous 
LD simulations and results relevant for thermonuclear flames. 

Sivashinsky (1977) and Michelson & Sivashinsky (1977) studied the non-linear regime 
of the LD instability analytically and numerically, using an integro-differential equation to 
describe the deformation of a discontinuous flame front. This was extended later by Gutman 
& Sivashinsky (1990), who concentrated on the regime where thermal conduction dominates 
over mass diffusion (this is the regime appropriate for astrophysical flames). They experi- 
ment with boundary conditions (both adiabatic and periodic) and discuss the influence of 
the domain size on the final cusp behavior. Their simulations show a multi-mode perturba- 
tion cusping, and slowly, the cusps merging until only a single one survives. In the widest 
domains, a superimposed cellular structure appears on the single remaining cusp. The mech- 
anism behind this cellular instability does not, at this moment, seem to be well understood 
(although see Kupervasser et al. 1996). A similar approach, using a more general equation 
for the front, the Frankel equation, considering spherically expanding flames, was performed 
by Blinnikov & Sasorov (1996). They observe a fractal structure for the front through a 
range of spatial scales. The unstable flame exhibits cell-splitting, owing to the expanding 
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spherical geometry. Helenbrook & Law (1999) presented simulations that treat the fuel and 
ash states as incompressible, and link them with jump conditions across the discontinuous 
front. In contrast to other methods, feedback between the flame and fluid flow is included 
in this model. They found that the flow dominates the LD instability in setting the scale for 
the wrinkling. This has important consequences for the flame in Type la supernova (SNe 
la), as the medium will be turbulent, and some evidence suggests that a large-scale dipole 
flow is setup during the ignition phase of the explosion (Kuhlen et al. 2003). In the real star, 
it may be likely that the fluid motions swamp any effects of the LD instability. Furthermore, 
the LD instability will be competing with the Rayleigh- Taylor instability (Bell et al. 2003b). 

Recently, front tracking methods treating a discontinuous flame coupled to full hydro- 
dynamics have been used for both terrestrial (Qian et al. 1998) and astrophysical flames 
(Ropke et al. 2003) to investigate LD unstable flames. The work of Ropke et al. (2003) used 
a level set method coupled to fully compressible PPM (Colella & Woodward 1984) to follow 
a discontinuous flame front, using flame speeds from Timmes & Woosley (1992). In their 
highest resolution run, they report seeing a cellular structure distorting the single cusp, per- 
haps arising through the same process as described in Gutman & Sivashinsky (1990). Ropke 
et al. (2003) find the flame speed rising to 1.3 x the laminar speed. 

The work described above all treat the flame as a discontinuous front — some model is 
required to describe the behavior of the flame on scales smaller than the grid resolution. The 
finite thickness of the flame is known to have some influence on the LD instability, such as 
setting the small scale cutoff for the growth, so simulations that resolve the thermal structure 
are important. Direct numerical simulations of the LD instability at various densities in 
the astrophysical context were performed by Niemeyer & Hillebrandt (1995), showing cusp 
formation, but it was unclear from their studies whether there was any acceleration over 
the laminar speed. Their lowest density flame (5 x 10 7 g cm" 3 ) showed significant growth 
beyond simple cusping, but it is likely that this was either numerical in origin, or due to the 
interaction of the transient perturbations left over from the initialization with the flame. 

LD unstable flames in laminar flows are unlikely alone to lead to the accelerations 
required of the SNe la models. An interesting question is whether the non-linear stabilization 
can break down, whether through some property intrinsic to the flames or through external 
forces. Kerstein (1996) discusses one such mechanism, active turbulent combustion, involving 
the interaction of the LD unstable flame with turbulence. It has been proposed that this 
process could occur in Type la supernovae (Niemeyer & Woosley 1997), since large amounts 
of the thermonuclear energy released goes into the expansion of the star, and it is this 
thermal expansion that drives the LD instability. Niemeyer & Woosley (1997) argue that 
this feedback between the flame and its own turbulence could lead to large accelerations of 
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the flame, although, at present, this has not been demonstrated in numerical simulations. 

We focus on DNS of the Landau-Darrieus instability at low densities (p ~ 2 x 10 7 g cm -3 
- 8 x 10 7 g cm" 3 ) using a low Mach number numerical formulation, as described in Bell 
et al. (2003a). Here we are able to extend to lower densities than previous studies, and 
also, through the use of adaptive mesh refinement, use finer resolution and larger domains. 
For the present study, we restrict ourselves to two-dimensional calculations. These direct 
numerical simulations are truly parameter free (i.e., no laminar flame speed is input), and 
could serve to test/calibrate different methods of representing unresolved flames on stellar- 
sized domains. Corresponding DNS calculations of the Rayleigh- Taylor instability will be 
presented in a second paper (Bell et al. 2003b). Three-dimensional results will be presented 
in a future paper. Additionally, we do not look to seed the flow with turbulence in the 
present study, but rather choose to focus solely on the pure LD instability. We provide a 
brief description of the numerics and parameters used for the present simulations in §2. In 
§3, we discuss the results, and finally, in §4 we conclude. 



2. NUMERICAL METHODS 

The evolution was carried out using a low Mach number formulation of the Euler equa- 
tions, as described in Bell et al. (2003a); Day & Bell (2000). The low Mach number for- 
mulation takes advantage of the fact that for slow moving flames, the pressure is constant 
across the flame front, to a high degree. Scaling the state variables and expanding them in 
powers of Mach number, keeping terms to 0(M 2 ) (Majda & Sethian 1985) yields a system 
of equations whose numerical stability condition is set by the fluid velocity alone, instead of 
the flow velocity plus sound speed. The pressure is decomposed into a dynamic and ther- 
modynamic component, the ratio of which is of 0(M 2 ), with only the dynamic component 
appearing in the momentum equation. The equation of state is enforced by requiring the 
pressure to be constant along particle paths, resulting in an elliptic constraint on the velocity 
field, completing the system of equations. For the flames we consider in the present paper, 
M < 10~ 4 , so this approximation is valid. This allows for timesteps that are 0(1/ M), larger 
than a corresponding compressible code. 

The low Mach number equations are solved via a second-order, fractional-step advec- 
tion/projection/reaction procedure. The advection step advances the fluid variables to the 
new timestep and finds a provisional velocity field that does not yet satisfy the elliptic con- 
straint. The new velocities are found by projecting this velocity field onto the space that 
discretely satisfies the elliptic constraint. Reactions are included via Strang-splitting. We 
refer the reader to Bell et al. (2003a) for the full details of the numerical method, as it is 
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unchanged for the present study. Laminar astrophysical flames computed with this code 
were compared to the benchmark calculations presented in Dursi et al. (2003), and found to 
be in good agreement. 

The calculations we present here differ from the fully compressible PPM (Colella & 
Woodward 1984) method used in Niemeyer & Hillebrandt (1995) in several fundamental 
ways. First, the low Mach number formulation yields a timestep constraint that depends 
only on the fluid velocity, not the fluid velocity and sound speed as in the compressible case. 
This means that we can take much larger timesteps, which for very slow moving flames, 
greatly reduces the accumulation of error. Additionally, the advection algorithm we employ 
is not dimensionally split, reducing the influence of the grid on the flame. Finally, the two 
dimensional domain is initialized by mapping a steady state flame in a perturbed fashion 
onto the grid, rather than perturbing an interface between fuel and ash and waiting for the 
flame to develop. This significantly reduces the transients and leaves the flow ahead of the 
flame unperturbed. 

An equation of state with degenerate/relativistic electrons and positrons, ideal gas ions, 
and Boltzmann radiation was used, as described in Timmes & Swesty (2000). The thermal 
conductivity is described in Timmes (2000). A single reaction, 12 C( 12 C,7) 24 Mg, is followed, 
using the unscreened rate from Caughlan & Fowler (1988). This accounts for the bulk of 
the energy generation in a C/O flame. Furthermore, at the densities we are considering, 
burning is expected to terminate at intermediate mass elements, so stopping at 24 Mg is 
reasonable. For this simple burning rate, we ignore the effects of the oxygen on the energy 
generation — it is only advected passively here. The temperatures reached behind the flame 
at these densities are ~ 4 x 10 9 K. At these temperatures, oxygen burning is several orders of 
magnitude slower than carbon, and we would not expect any effects from it on the timescales 
that we follow the flames. The properties of the flames we consider here are listed in Table 1. 
We will consider the 4 x 10 7 g cm" 3 flame to be our reference flame, upon which we will 
explore the effects of resolution and domain size. 

Motivated by results indicating that multi-mode perturbations eventually merge into a 
single dominant cusp (see for example Gutman & Sivashinsky 1990; Bell et al. 2003a), we 
concentrate mainly on single mode perturbations. In all cases, the flame is initialized on the 
grid by mapping a steady-state one-dimensional laminar flame across the grid, shifting it 
vertically according to the perturbation. The boundary conditions in all cases are periodic 
transverse to the flame propagation direction, inflow ahead of the flame, and outflow behind 
the flame. The flame is oriented such that it is propagating downward in our domain. The 
inflow conditions correspond to the fuel conditions, with the inflow velocity chosen to match 
the laminar flame velocity. Thus, a flame moving at the laminar speed will remain stationary 
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in the computational domain. Any acceleration resulting from the flame instability will lead 
to motion across the grid. The resolution used is chosen to put about 5-10 zones inside the 
thermal width, If, calculated as 

j A ash A fuel / „ \ 

f = max(VT) ' 1 j 

where T asn and Tf ue i are the ash and fuel temperatures respectively. We note that this flame 
width measure is ~ 2 times narrower than that used in Timmes & Woosley (1992) (see also 
Dursi et al. 2003). In §3, we look at the sensitivity of the results to the resolution. 



3. RESULTS 

3.1. Multi-mode perturbations 

As a motivation for the single mode studies, we look at a single case of a multi-mode 
perturbation of a 4 x 10 7 g cm" 3 flame in a 20.48 cm wide domain. The flame front was 
perturbed by shifting the zero point of the initial steady state flame by a sinusoid, with 30 
different frequencies, whose phases and amplitudes were chosen randomly. The maximum 
amplitude of the perturbation is one flame thickness. Figure 1 shows the y- velocity at several 
instances in time. We see that immediately, the higher frequency perturbations are washed 
out, and only the long wavelength modes slowly begin to grow. Several cusps begin to form, 
but quickly merge, leaving two dominant cusps behind. These remaining cusps gradually 
move toward each other (through the periodic boundary), eventually merging, leaving a 
single dominant cusp behind. The timescale for the cusps to merge is long, ~ 2 ms, so we 
show only one case here. Bell et al. (2003a) shows another example, for a pure carbon flame 
at a higher density, following the evolution for several cusp mergers. 

The flame speed can be measured by computing the carbon destruction rate on the grid, 
taking into account the flow of material through the boundaries: 

F(() = - wtir (4) 

where Q is the spatial domain of the burning region, W is the width of inflow face, (pX c ) m 
is the inflow carbon mass fraction, and puj c is the rate of consumption of carbon due to 
nuclear burning. The rate of change of this quantity in a time interval [Ti,T2] is a measure 
of effective flame velocity: 

eff_ {T 2 - Tl )w{ P x c r {) 
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Figure 2 shows the flame speed as a function of time for the multimode calculation. After 
an initial, relatively constant velocity period, in which the initial perturbations die out, the 
flame begins to accelerate rapidly. The velocity peaks shortly after the leftmost two initial 
cusps merge, leaving behind two well defined cusps. A second peak occurs at the merging of 
these two cusps. 

The merging of multiple cusps into a single dominant one has been observed in other 
planar flame simulations (see for example Gutman & Sivashinsky 1990; Helenbrook & Law 
1999), and it seems to be a general feature of the LD instability, with the exception of 
spherically expanding flames (Blinnikov & Sasorov 1996). We note however that we see no 
evidence for a superposed cellular structure on the remaining cusp, but this may be due 
to domain size restrictions (Gutman & Sivashinsky 1990). The final cusp that emerges is 
robust and shows no signs of breaking down. 

3.2. Single-mode perturbations 

The main focus of the present study is to look at the acceleration that a flame undergoing 
the LD instability will experience, in the absence of any external forcing. In particular, we 
wish to determine whether there is any breakdown of the cusp configuration. Since the 
general consensus is that multimode perturbations in planar domains tend to merge into a 
single, well defined cusp, we can seed a single mode to study the behavior of this steady 
state configuration. 

Table 2 list the parameters used for the single-mode calculations. In addition to a run 
of densities, we considered the effect of the box width and resolution on the instability for 
one choice of density (4 x 10 7 g cm -3 ). The reference calculation at this density is a 10.24 cm 
wide x 20.48 cm high domain, with a resolution of 6 points in the flame width, If. Numerical 
experiments have shown that this resolution is enough to accurately follow the flame front 
(see Bell et al. 2003a). The critical wavelength for the growth of 1/2 carbon, 1/2 oxygen 
astrophysical flames over a range of densities was computed in Dursi et al. (2003), who, 
when converted to the flame thickness definition used here, found ~ 80 flame thicknesses 
to be the small scale cutoff. If this result holds for the densities we consider here, then for 
domains narrower than 80 flame widths, we do not expect the LD instability to be able to 
grow. We note that our reference domain width is ~ 174 flame widths wide, which should 
be well above the small scale cutoff for the LD growth. We will find this small scale cutoff 
directly by varying our domain width, computing the Markstein number, and finding the 
zeros of the LD dispersion relation, Equation (1). A single sine wave perturbation was used 
as an offset when mapping the steady state laminar flame profile onto the domain. 
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Figure 3 shows the x- and y- velocities for the reference 4 x 10 7 g cm -3 LD flame once the 
cusp has become well defined. After this point, the flame front retains this shape. As with 
the cusps in the multimode simulation, we see that the vertical velocity peaks strongly right 
behind the cusp, and then drops below the ambient post flame velocity as one moves further 
away from the flame front. One thing to note from these velocity figures is how smooth 
the flow is ahead of and behind the flame. There is virtually no noise to disturb the pure 
LD instability. Figure 4 illustrates this by showing two slices of the ^-velocity — through 
the center of the cusp and the center of the trough — when the cusp is at its maximum. 
This behavior is a natural consequence of the local curvature and its effect on the velocity 
streamlines. The change in the peak carbon destruction rate as a function of curvature is 
shown in Figure 5. The curvature was computed by finding the position of the front and 
taking the second derivative. We see that the larger the magnitude of the curvature, the 
greater the carbon destruction rate. At zero curvature, the carbon destruction rate is the 
laminar value. The scatter in the plot reflects the difficultly in defining the curvature for a 
discretely defined interface. 

Figure 6 shows the flame speed as a function of time for the 4 x 10 7 g cm~ 3 C/O flame in 
several different sized domains (2.56 cm, 5.12 cm, 10.24 cm, and 20.48 cm). The flame speed 
was computed by looking at the carbon consumption rate as in Equation (5). As the box 
width was varied, we choose to make 5y/W constant, where Sy is the perturbation amplitude, 
and W is the box width. All other parameters (i.e., box height, boundary conditions, etc.) 
were held constant across the different simulations. The velocity of the flame for the two 
widest domains increases quickly and reaches a peak, where it levels off, at a value about 
1.6% higher than the laminar speed. The slight decrease in the flame speed at late times for 
the 10.24 cm run is likely a result of the finite vertical extent of the domain. The 5.12 cm 
domain run evolves much more slowly. This domain is ~ 87 flame widths wide, right at the 
small scale cutoff for the LD instability. We believe that the slow growth that we see is a 
manifestation of this cutoff. The narrowest domain (2.56 cm) does not grow at all. It is 
not known what accounts for the small difference in the asymptotic flame speeds for the two 
widest domains. 

We can compare the time it takes for the flame speed to saturate to the theoretical 
prediction of the LD growth rate, Equation (1). The amplitude of the cusp is computed by 
laterally averaging the carbon mass fraction, and finding the positions where it exceeds 0.05 
and 0.45. Figure 7 shows the height of the cusp as a function of time for the 10.24 cm 
and 20.48 cm wide runs. Looking at these plots, there is clearly a regime of linear growth 
followed by a transition to the nonlinear regime where the cusping halts the growth, and the 
amplitude reaches a steady state value. The solid line is a fit to 

V = Ae ui , (6) 
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for the interval of 1 x 10~ 4 to 4 x 10 -4 s (the linear regime), which gives uj = 3968 and 
3046 respectively — the narrower domain grows faster. We ignored the very first part of the 
growth to avoid any transients. If we were to ignore the effects of curvature and the flame 
thickness (i.e., set Ma = in Equation [1]), we would get predictions of u — 8200 and 4100 
for these domains. Thus, the finite thickness of the flame slows down the growth of the LD 
instability, demonstrating the effects of the curvature on the flame speed. 

The finite thickness imposes a small scale cutoff to the LD growth, which we can now 
compute. Using the flame data from Table 1, we can solve for the Markstein number, and 
find Ma = —2.45 for both of these domain widths. A negative value is to be expected, as it 
indicates that a positively curved flame burns more slowly, which is expected for our flames 
where thermal diffusion dominates mass diffusion (Dursi et al. 2003). With this value of the 
Markstein number, we can plot the growth rate as a function of wavenumber, Figure 8. We 
see that the LD stops growing for dimensionless wavenumbers larger than 0.07, or domain 
widths smaller than 5.28 cm. This is consistent with behavior of the 2.56 cm and 5.12 cm 
runs we showed above. For this reason, we do not perform fits to these narrowest domain 
runs. The maximum growth should be seen for dimensionless wavenumbers of 0.035, or a 
domain with of 10.6 cm. Looking back at Figure 6, we see that the 10.24 cm run (our closest 
domain width to this maximum) reached the peak velocity quickest and had the largest 
growth rate (as determined by the fits to the cusp amplitude growth), consistent with the 
dispersion relation. This also explains why the multimode perturbation presented above went 
through an intermediate phase consisting of two cusps. There, the domain was 20.48 cm 
wide, and perturbations were made on many length scales. The dispersion relation predicts 
that the 10.6 cm perturbations will grow fastest, overwhelming any smaller perturbations. 
This means that the mode that we seeded that put two wavelengths in the box would grow 
the fastest, resulting in the intermediate, moderately long lived state with two cusps. 

Figure 9 shows the flame speed vs. time for the 10.24 cm wide, 4 x 10 7 g cm -3 run 
at three different resolutions. The flame speed is normalized to the laminar speed at that 
resolution to account for the small resolution dependence on this speed. Table 3 lists the 
laminar flame speeds for the different resolutions. At the highest resolutions, the laminar 
flame speeds differ by 0.13%, and the normalized velocity vs. time curves lie on top of one 
another, demonstrating convergence. The coarsest resolution laminar flame speed differs by 
1.1%, and the evolution with time is clearly not converged. All the simulations presented 
in this paper are performed at this middle resolution or higher. This resolution agrees with 
that found in the corresponding resolution study in Bell et al. (2003a). 

Figures 10 and 11 show the flame speed as a function of time for the 2 x 10 g cm -3 and 
8 x 10 7 g cm -3 runs respectively. The evolution at these densities proceeded in much the same 



-10- 



way as the 4 x 10 7 g cm -3 run shown in Figure 3, so we do not show the velocity fields. Only 
one domain width was used for these calculations, with the size picked close to the wavelength 
that maximizes the growth rate. The change in density (and expansion parameter across 
the flame) changes the speedup we get. The 2 x 10 7 g cm -3 flame accelerates by 2.4% and 
the 8 x 10 7 g cm -3 flame by 1.1%. Together with the results for the 4 x 10 7 g cm -3 flame, 
we see that the acceleration is larger for the lower density/larger expansion ratio flames. 
We can compute the Markstein numbers for these flames in the same manner as above and 
find —3.42 and —2.02 for the 2 x 10 7 g cm -3 and 8 x 10 7 g cm -3 flames respectively (see 
Figures 12 and 13). Again, the fit of an exponential to the linear portion of the growth is 
excellent. The dispersion relations for these flames are also plotted in Figure 8. To fit them 
all on the same plot, the curves are scaled by the flame speed and width. The magnitude of 
the Markstein number increases with decreasing density and increasing expansion parameter. 
This is not unexpected, and confirms the trend seen in Dursi et al. (2003) for similar flames. 
This suggests that a lower density flame is harder to bend. 

As a final measure of the effects of curvature on the flame speed, we can compare the 
growth in the flame surface area to the growth in flame speed. In the simplest model, these 
two are directly related in a purely geometrical fashion, and the flame speed would be simply 



We can measure the length of the flame surface by counting the number of zone edges where 
the carbon mass fraction passes through 0.25. This measure agrees well with computing 
the length of a contour as computed with the commercial IDL package, after normalization. 
Since the flame surface is always sharp, this is a well-defined procedure. Figure 14 shows the 
normalized flame length as a function of time for the three different densities (the 10.24 cm 
reference run was used for the 4 x 10 7 g cm -3 case). In all these cases, the increase in the 
flame length is several times larger than the increase in the flame speed, showing a strong 
deviation from the linear scaling of Equation (7) — for example, the 4 x 10 7 g cm" 3 flame 
increases in speed by about 1.6%, but the area increase is 12%. This deviation is due to the 
effects of the localized curvature. We note that this deviation is larger in magnitude than 
that found in the compressible, flame- model calculations of Ropke et al. (2003). 



v(t) = v 



A(t) 



(7) 



A 



4. 



CONCLUSIONS 



We presented direct numerical simulations of the Landau-Darrieus instability at low 
densities in conditions relevant for Type la supernova explosions. We employed a low Mach 
number hydrodynamics method — new in the application to astrophysics — to advance the 
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flows and enable us to follow the development of the LD instability from the linear regime 
through the nonlinear regime. Resolving the flame ensures that we account for the full 
coupling between the flame and the flow and the effects of curvature. The consistency of our 
results with the theoretical predictions validate this approximation for astrophysical flames. 
Fully compressible versions of the same simulations would not be possible, or prohibitively 
expensive and error prone. 

Acceleration of the flame was seen, due to the increase of flame surface, until the non- 
linear cusp formation sets in, at which point the velocity stabilized a few percent higher 
than the laminar speed. The maximum increase in flame speed observed was ~ 2% — this 
is in contrast to the results presented in Ropke et al. (2003), where accelerations of ~ 30% 
were observed. There are several possible reasons to explain these differences. Most of the 
numerical results for the LD instability treat the flame as infinitely thin. Here we've included 
the effects of the finite thickness of the flame, and therefore, implicitly included curvature ef- 
fects. At the densities we consider, the burning is expected to terminate at the intermediate 
mass elements, so we burn to magnesium whereas they burn to nickel. Furthermore we use 
half carbon/half oxygen fuel instead of pure carbon. This will affect the expansion factor 
across the flame. Finally, fully understanding the role of the numerics requires a detailed 
comparison of the different methods, using the same input physics, and could potentially be 
a research project in the future. 

The speedup increases with decreasing density (or increasing expansion parameter), 
as expected. With time, it seems as if the velocities asymptote to the same value, for a 
given density. When a multimode perturbation is considered, it is observed that the cusps 
eventually merge into a single cusp. 

We were able to perform fits of the growth of the cusp amplitude with time to compute 
the growth rates. Using the theoretical prediction (Equation [1]), we were able to compute 
the Markstein number for these flames, finding values < —2.0, and growing in magnitude 
with decreasing density. These values and the trend with density are consistent with the 
results of Dursi et al. (2003). For the 4 x 10 7 g cm -3 density run, the large number of domain 
widths allowed us to compare to the predictions of the fastest growing mode and small scale 
cutoff from the theoretical dispersion relation, and we found strong agreement. Furthermore, 
this demonstration of the small scale cutoff and the deviation from the infinitely thin flame 
(Ma = 0) estimates of the growth rate, as well as the difference in the growth of the flame 
speed and area emphasize the importance of curvature in small-scale flame dynamics. As 
we zoom out from these fully resolved flames, several steps will be needed at intermediate 
scales between the fully resolved flame, where curvature effects appear important, and the 
full star, with validation of the flame model on each scale through direct comparison to the 
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immediately smaller scale simulations. The present direct numerical simulations provide the 
anchor for this subgrid model validation. These curvature effects will be explored more fully 
in paper II, focusing on the reactive Rayleigh- Taylor instability. 

We saw no breakdown of the nonlinear stabilization mechanism as hinted to in the 
calculations presented by Niemeyer & Hillebrandt (1995). Furthermore, none of the flames 
we studied showed a superposed cellular structure, although, there is some evidence (Gutman 
& Sivashinsky 1990) that this phenomena is present only on very large scales, and Blinnikov 
& Sasorov (1996) admit that it may actually be seeded by numerical noise. The widest 
domain that we were able to simulate was ~ 350 flame widths. Helenbrook & Law (1999) 
only saw this instability when they increased their domain size to 20 times the wavelength 
where the LD growth rate is the maximum. This would translate into a domain size of 
230 cm for the 4 x 10 7 g cm -3 flame, which would require a 23,000 zone wide simulation if 
we are to resolve the flame. The growth rate for a domain that wide would be very small, 
requiring a much longer simulation. This is very expensive, and furthermore, it is not known 
if this factor of 20 is universal. We are not aware of any resolved calculations of the LD 
instability that have shown this secondary instability. 

We note that all of our calculations considered planar geometries only. On the scales 
where the spherical geometry of the star is likely to make a difference, the flame will be 
disturbed by the larger scale motions in the star, other instabilities as well (notably the 
Rayleigh- Taylor instability, presented in Bell et al. 2003b) and subjected to the interaction 
with flame-generated turbulence. All of these processes will dominate the pure LD instability, 
considered in the present paper. 

The authors thank L. J. Dursi for useful comments on the manuscript and F. X. Timmes 
for making his equation of state and conductivity routines available on-line. Support for this 
work was provided by the DOE grant number DE-FC02-01ER41176 to the Supernova Science 
Center/UCSC and the Applied Mathematics Program of the DOE Office of Mathematics, 
Information, and Computational Sciences under the U.S. Department of Energy under con- 
tract No. DE-AC03-76SF00098. The calculations presented here were performed on the 
IBM Power4 (cheetah) at ORNL, sponsored by the Mathematical, Information, and Com- 
putational Sciences Division; Office of Advanced Scientific Computing Research; U.S. DOE, 
under Contract No. DE-AC05-00OR22725 with UT-Battelle, LLC and the UCSC UpsAnd 
cluster supported by an NSF MRI grant AST-0079757. 



13 



REFERENCES 

Bell, J. B., Day, M. S., Rendleman, C. A., Woosley, S. E., & Zingale, M. A. 
2003a, Journal of Computational Physics, to appear, preprint available at: 
http: / / seesar.lbl.gov/CCSE/Publications/car/LMN_nucl_flames.pdf 

— . 2003b, ApJ, in preparation 

Blinnikov, S. L, & Sasorov, P. V. 1996, Phys. Rev. E, 53, 4827 

Caughlan, G. R., & Fowler, W. A. 1988, Atomic Data and Nuclear Data Tables, 40, 283, see 
also http:/ /www. phy.ornl.gov/astrophysics/data/cf88/index. html 

Clanet, C, & Searby, G. 1998, Physical Review Letters, 80, 3867 

Colella, P., & Woodward, P. R. 1984, J. Comp. Phys., 54, 174 

Darrieus, G. 1938, unpublished work presented at La Technique Moderne 

Day, M. S., & Bell, J. B. 2000, Combust. Theory Modelling, 4, 535 

Dursi, L. J., Zingale, M., Calder, A. C, Fryxell, B., Timmes, F. X., Vladimirova, N., Rosner, 
R., Caceres, A., Lamb, D. Q., Olson, K., Ricker, P. M., Riley, K., Siegel, A., & Truran, 
J. W. 2003, ApJ, 595, 955 

Gamezo, V. N., Khokhlov, A. M., Oran, E. S., Chtchelkanova, A. Y., & Rosenberg, R. O. 
2003, Science, 299, 77 

Gutman, S., & Sivashinsky, G. I. 1990, Physica D, 43, 129 

Helenbrook, B. T., & Law, C. K. 1999, Comb. Flame, 117, 155 

Hillebrandt, W., & Niemeyer, J. C. 2000, Annu. Rev. Astron. Astrophys, 38, 191 

Kerstein, A. R. 1996, Combustion Sci. Technol., 118, 189 

Kuhlen, M., Woosley, S. E., & Glatzmaier, G. 2003, Astrophysical Journal, in preparation 

Kupervasser, O., Olami, Z., & Procaccia, I. 1996, Physical Review Letters, 76, 146 

Landau, L. D. 1944, Acta Physicochimica, USSR, 19, 77, reprinted in P. Pelce, Dynamics of 
Curved Fronts, Academic Press, Berkeley, CA, 1988 

Majda, A., & Sethian, J. A. 1985, Combust. Sci. Technol., 42, 185 



-14- 



Michelson, D. M., & Sivashinsky, G. I. 1977, Acta Astronautica, 4, 1207 

Niemeyer, J. C, & Hillebrandt, W. 1995, ApJ, 452, 779 

Niemeyer, J. C, & Woosley, S. E. 1997, ApJ, 475, 740 

Qian, J., Tryggvason, G., & Law, C. K. 1998, JCP, 144, 52 

R6pke, F. K., Niemeyer, J. C, & Hillebrandt, W. 2003, ApJ, 588, 952 

Reinecke, M., Hillebrandt, W., & Niemeyer, J. C. 2002, A&A, 391, 1167 

Sivashinsky, G. I. 1977, Acta Astronautica, 4, 1177 

Timmes, F. X. 2000, ApJ, 528, 913 

Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 

Timmes, F. X., & Woosley, S. E. 1992, ApJ, 396, 649 

Zeldovich, Y. B., Barenblatt, G. I., Librovich, V. B., & Markhviladze, G. M. 1985, The 
Mathematical Theory of Combustion and Explosions (New York: Consultants Bu- 
reau) 



This preprint was prepared with the AAS IATgX macros v5.2. 



- 15 - 



Table 1: Flame properties for 0.5 12 C/0.5 16 flames. 



p 


a a 


^laminar 


h 


(g cm -3 ) 




(cm s" 1 ) 


(cm) 


2 x 10 7 


1.69 


1.50 x 10 4 


0.26 


4 x 10 7 


1.52 


6.08 x 10 4 


0.059 


8 x 10 7 


1.39 


2.04 x 10 5 


0.011 



a a = Pfuel/Pash 



Table 2: Parameters for single mode Landau-Darrieus simulations. 



p 


description 


W x H 


fine g 


;rid 


W/lf 


l f /Ax 


(g cm -3 ) 




(cm) 










2 x 10 7 




51.20 x 


102.4 


1536 x 


3072 


199.2 


7.7 


4 x 10 7 


very narrow 


2.56 x 


20.48 


256 x 


2048 


43.4 


5.9 




narrow 


5.12 x 


20.48 


512 x 


2048 


86.8 


5.9 




regular 


10.24 x 


20.48 


1024 x 


2048 


173.6 


5.9 




wide 


20.48 x 


20.48 


2048 x 


2048 


347.2 


5.9 




low-resolution 


10.24 x 


20.48 


512 x 


2048 


173.6 


3.0 




high-resolution 


10.24 x 


20.48 


2048 x 


4096 


173.6 


11.8 


8 x 10 7 




2.56 x 


3.84 


2048 x 


3072 


232.7 


8.8 



Table 3: 4 x 10 7 g cm 3 laminar flame speed vs. resolution 



If /Ax 


^laminar 




(cm s _1 ) 


3.0 


6.013 x 10 4 


5.9 


6.074 x 10 4 


11.8 


6.082 x 10 4 
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Fig. 1. — y- velocity for a 4 x 10 7 g cm -3 C/O flame undergoing a Landau-Darrieus instability. 
In this simulation, 30 different modes were perturbed. Each pane is 3x 10~ 4 s apart, spanning 
s to 2.1 x 10~ 3 s. The initial perturbation rapidly evolves to only a few growing modes, and 
quickly, two cusps form that dominate the flow. Slowly, these two cusps move toward each 
other and merge through the periodic boundary, leaving a single, well-defined cusp behind. 
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Fig. 2. — Flame speed as a function of time for the 4 x 10 7 g cm 3 C/O multimode Landau- 
Darrieus instability, shown in Figure 1. 
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Fig. 3. — x-velocity (left) and ^-velocity (right) for the high-resolution, 10.24 cm wide 
4 x 10 7 g cm -3 flame perturbed with a single mode after 5 x 10 -4 s. 
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Fig. 4. — y-velocity for the 10.24 cm wide, 4 x 10 7 g cm" 3 flame perturbed when the flame 
velocity is maximum (6.29 x 10~ 4 s) at the peak of the cusp (solid line) and the middle of the 
trough (dotted line). This demonstrates that the curvature leads to a change in the flame 
speed at different points along the flame front. 
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Fig. 5. — Carbon destruction rate vs. curvature (rf 2 |/interface/^a ;2 ) for the 10.24 cm wide, 
4 x 10 7 g cm" 3 flame. When there is zero curvature, the carbon destruction rate is the 
laminar value. 
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Fig. 6. — Flame speed vs. time for the 4 x 10 7 g cm -3 LD unstable flame in a 20.48 cm 
domain (dot) 10.24 cm domain (dash), 5.12 cm domain (solid), and 2.56 cm domain (dot- 
dash). The slow growth of the 5.12 cm wide run reflects the influence of this width being 
near the small scale cutoff for the growth of the LD instability. The 2.56 cm run shows no 
signs of growing, as expected, since that domain width is well below the small scale cutoff. 
The 10.24 cm wide domain rises the fastest, as expected, since this is closest to the maximum 
u{k) of Equation [1] for this flame. 
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Fig. 7. — Height of the cusp (natural log scale) as a function of time (symbols) and a fit on 
an exponential to the linear part (solid line) for the 10.24 cm wide domain (left) and the 
20.48 cm wide domain. 
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0.005 




Fig. 8. — Predicted non-dimensional growth rate (u lf/Uj) of the LD instability as a function 
of dimensionless wavenumber (k If) from Equation [1] for the three different densities, 2 x 
10 7 g cm -3 (dashed), 4 x 10 7 g cm -3 (solid), and 8 x 10 7 g cm -3 (dotted), based on the 
measured Markstein numbers. 
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0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 
time(s) 

Fig. 9. — Normalized flame speed, Ves(t)/V e s(t = 0), vs. time for the 4 x 10 7 g cm -3 
LD unstable in a 10.24 cm wide domain at 3 different resolutions: 3.0 zones/// (solid), 
5.9 zones/// (dash), and 11.8 zones/// (dot). The two highest resolution runs seem to have 
converged. The main calculations presented in this paper were computed at the 5.9 zones/// 
resolution. The normalization is required to correct for the small differences in the laminar 
flame speed with resolution. 



-25 - 




Fig. 10. — Flame speed as a function of time for the 2 x 10 7 g cm 3 LD unstable flame. 
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Fig. 11. — Flame speed as a function of time for the 8 x 10 7 g cm 3 LD unstable flame. 
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Fig. 12. — Cusp amplitude (natural log scale) as a function of time (symbols) and a fit on 
an exponential to the linear part (solid line) for the 2 x 10 7 g cm -3 flame, in the interval 
[1CT 3 s, 4 x 10~ 3 s]. The growth rate of the fit is u = 244 s _1 . 
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Fig. 13. — Cusp amplitude (natural log scale) as a function of time (symbols) and a fit on 
an exponential to the linear part (solid line) for the 8 x 10 7 g cm -3 flame, in the interval 
[5 x 10~ 6 s, 2 x 10" 5 s]. The growth rate of the fit is u = 522000 s -1 . 
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2 x 1 7 g cm 3 4 x 1 7 g cm 3 8 x 10' g cm 3 




0.000 0.004 0.008 0.012 0.016 0.0000 0.0002 0.0004 0.0006 0.0008 2x10-5 4x10-5 6x10-5 

time(s) time(s) time(s) 

Fig. 14. — Normalized flame length, L(t)/L(t = 0), as a function of time for the three 
different densities, 2 x 10 7 g cm -3 (left), 4 x 10 7 g cm -3 (center), and 8 x 10 7 g cm -3 (right). 
The increase in flame length is several times larger than the corresponding increase in the 
flame speed. 



